function stokes_test

addpath ../

p = params;

p.Ma = 1;

global N Nz k dz

k_min = 1;

kk = 3*k_min;

N = 2^6;
Nz = 50;

h = 2 * pi / N;
dz = 1 / (Nz - 1);

k = k_min * [0:N/2 - 1 0 -N/2+1:-1];

x = h * (1:N) / k_min;
z = linspace(0, 1, Nz);

c_top = sin(3*k_min * x);

c_hat = fft(c_top);
h = 1;
for j = 1:Nz
    w_hat = -0.4e1 * p.Ma * exp(k) .* c_hat .* (((-0.1e1 / 0.2e1 + k * z(j)) .* (exp(0.2e1 * k) / 0.2e1 + k + 0.1e1 / 0.2e1) .* exp(0.2e1 * k * z(j)) + ((-z(j) / 0.2e1 + k * z(j) + 0.1e1 / 0.2e1) .* exp(0.2e1 * k) + 0.1e1 / 0.2e1 - z(j) / 0.2e1) .* k) .* exp(-k * z(j)) - exp(k * z(j)) .* (-0.1e1 / 0.2e1 + (-0.1e1 / 0.2e1 + k) .* exp(0.2e1 * k)) / 0.2e1) ./ (0.4e1 * exp(0.4e1 * k) + 0.8e1 * exp(0.2e1 * k) + 0.16e2 .* k .^ 2 .* exp(0.2e1 * k) + 0.4e1);
%     w_hat = -0.2e1 * exp(k * h) .* c_hat .* (((-(p.Ma + 2 * p.delta) .* (k * z(j) - 0.1e1 / 0.2e1) .* exp(0.2e1 * k * z(j)) - 0.2e1 * h * z(j) * (p.Ma + 2 * p.delta) * k .^ 2 + ((-0.2e1 * z(j) - 0.2e1 * h) * p.delta - p.Ma * (-z(j) + h)) .* k - (2 * p.delta)) .* exp(0.2e1 * k * h) + 0.2e1 * (k * z(j) - 0.1e1 / 0.2e1) .* (h * (p.Ma + 2 * p.delta) .* k + p.Ma / 0.2e1 - p.delta) .* exp(0.2e1 * k * z(j)) + (p.Ma + 2 * p.delta) * (-z(j) + h) .* k - (2 * p.delta)) .* exp(-k * z(j)) + ((h * (p.Ma + 2 * p.delta) * k - p.Ma / 0.2e1 + p.delta) .* exp(0.2e1 * k * h) + p.delta + p.Ma / 0.2e1) .* exp(k * z(j))) ./ (-0.4e1 * exp(0.4e1 * k * h) + 0.16e2 * k .* h .* exp(0.2e1 * k * h) + 0.4e1);
    w(j,:) = real(ifft(w_hat));
    u_hat = (-4*1i) * p.Ma * (((-0.1e1 / 0.2e1 + k * z(j)) .* (exp(0.2e1 * k) / 0.2e1 + k + 0.1e1 / 0.2e1) .* exp(0.2e1 * k * z(j)) + (-0.1e1 / 0.2e1 - k .^ 2 * z(j) + (0.1e1 / 0.2e1 + z(j) / 0.2e1) .* k) .* exp(0.2e1 * k) - 0.1e1 / 0.2e1 + (-0.1e1 / 0.2e1 + z(j) / 0.2e1) .* k) .* exp(-k * z(j)) - ((k - 0.3e1 / 0.2e1) .* exp(0.2e1 * k) - 0.2e1 * k - 0.3e1 / 0.2e1) .* exp(k * z(j)) / 0.2e1) .* exp(k) .* c_hat ./ (0.4e1 * exp(0.4e1 * k) + 0.8e1 * exp(0.2e1 * k) + 0.16e2 * k .^ 2 .* exp(0.2e1 * k) + 0.4e1);
%     u_hat = (-2*1i) .* exp(k * h) .* c_hat .* (((-(p.Ma + 2 * p.delta) .* (k * z(j) - 0.1e1 / 0.2e1) .* exp(0.2e1 * k * z(j)) + 0.2e1 * h * z(j) * (p.Ma + 2 * p.delta) * k .^ 2 + ((-h - z(j)) * p.Ma - 0.2e1 * (-z(j) + h) * p.delta) .* k + p.Ma) .* exp(0.2e1 * k * h) + 0.2e1 * (k * z(j) - 0.1e1 / 0.2e1) .* (h * (p.Ma + 2 * p.delta) .* k + p.Ma / 0.2e1 - p.delta) .* exp(0.2e1 * k * z(j)) - (p.Ma + 2 * p.delta) * (-z(j) + h) * k - p.Ma) .* exp(-k * z(j)) + exp(k * z(j)) .* ((h * (p.Ma + 2 * p.delta) .* k - 0.3e1 / 0.2e1 * p.Ma - p.delta) .* exp(0.2e1 * k * h) + 0.2e1 * h * (p.Ma + 2 * p.delta) * k + 0.3e1 / 0.2e1 * p.Ma - p.delta)) ./ (-0.4e1 * exp(0.4e1 * k * h) + 0.16e2 * k .* h .* exp(0.2e1 * k * h) + 0.4e1);
    u(j,:) = real(ifft(u_hat));
end

quiver(x, z, u, w);
norm(w(Nz,:) - p.delta * c_top)